/*
 *  viscoelastic.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 1/04/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "viscoelastic.h"

Viscoelastic::Viscoelastic()
:	Material(),
	vTau(),
	vBeta(),
	nTau(),
	sumbeta(0.0),
	rheological(false),
	matHyper()
{}

// ***

Viscoelastic::Viscoelastic(ifstream &in, const string& na, const string & t, Material* ma)
:	Material(na,t),
	vTau(),
	vBeta(),
	vMu(),
	vEta(),
	nTau(),
	sumbeta(0.0),
	integrated(true),
	rheological(false),
    Wrhe(0.0),
	matHyper(ma)
{
	int i;
	string keyword;
	string line;
	vector<string>	tauValues;
	stringstream stnum;
	double value;
	
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("RHEOLOGICAL"))
        {
			rheological = true;
        }
        //in >> keyword;
        if(!keyword.compare("PREINTEGRATED"))
        {
        	integrated = false;
        }
        //in >> keyword
        if(!keyword.compare("TAU"))
		{
        	getline(in,line);
            istringstream iss(line);
            while (iss >> keyword )
            {
            	tauValues.push_back(keyword);
            }
            nTau = tauValues.size();
			
            for(i=0; i<nTau; ++i)	// Starts from 1 because the first element is the keyword V0
            {
            	stnum << tauValues[i];
                stnum >> value;
                vTau.push_back(value);
                stnum.clear();
            }
		}
        if(!keyword.compare("BETA"))
        {
        	for (i=0; i<nTau; ++i)
        	{
        		in >> value;
        		vBeta.push_back(value);
        	}
        }
    }
    
    hypertype = matHyper->hypertype;

    vEta.setZero(nTau);
    vMu.setZero(nTau);
        
    if (integrated)
    {
    	sumbeta = 0.0;
    	for(i=0;i<nTau;++i)
        {
    		sumbeta += vBeta[i];
    		vMu[i] = vBeta[i] * matHyper->getE() / 2.0;
    		vEta[i] = 2.0 * vMu[i] * vTau[i];
        }
    }
    if (rheological)
    {
    	Idv.setZero(6);
        Idv.head(3) <<  1.0,
                        1.0,
                        1.0;
            
        Id4m = MatrixXd::Identity(6,6);
        Id4m.block(3,3,3,3) *= 0.5;
     }

    E.setZero(6);
    PMat.setZero(6,6);
    PtMat.setZero(6,6);
    tildePMat.setZero(6,6);
    tildeGamoCi.setZero(6,6);
    CiotildeGam.setZero(6,6);

    density = matHyper->getDensity();
}

// ***

void Viscoelastic::calcW( const VectorXd & C , const VectorXd & Gamma )
{
    matHyper->calcW(C);		// Here the invariants are computed for C
    W = matHyper->getW();
    
    switch (rheological)
    {
        case false:
            if(!hypertype.compare("coupled"))
            {
            	matHyper->calcS(C);
                for(int i=0 ; i<nTau ; ++i)
                    W += vMu[i] * dotVoigt(Gamma.segment<6>(6*i), Gamma.segment<6>(6*i)) - vBeta[i] * dotVoigt(matHyper->getS(), Gamma.segment<6>(6*i)) + vBeta[i] * matHyper->getW();
            }
            else
            {
            	matHyper->calcSiso(C);
                for(int i=0 ; i<nTau ; ++i)
                    W += vMu[i] * dotVoigt(Gamma.segment<6>(6*i), Gamma.segment<6>(6*i)) - vBeta[i] * dotVoigt(matHyper->getSiso(), Gamma.segment<6>(6*i)) + vBeta[i] * matHyper->getWiso();
                
            }
            break;
                
        case true:

            for(int i=0 ; i<nTau ; ++i)
            {
                calcWrhe(i,C);
            	calcSrhe(i,C);
                W += vMu[i] * dotVoigt(Gamma.segment<6>(6*i), Gamma.segment<6>(6*i)) - dotVoigt(Srhe, Gamma.segment<6>(6*i)) + Wrhe;
            }
            break;
    }
}

// ***

void Viscoelastic::calcS( const VectorXd & C, const VectorXd & Gamma )
{
    matHyper->calcS(C);
    S = matHyper->getS();

    if(type.compare("svk")) matHyper->calcCt(C); //Means that hyperelastic material is not svk.
                
    switch (rheological)
    {
        case false:
        
            if(!hypertype.compare("coupled"))
            {
                for(int i=0 ; i<nTau ; ++i)
                {
                	//S += vBeta[i] * ( matHyper->getS() - matHyper->getCt() * makeVoigt(Gamma.segment<6>(6*i)) );
                    S += vBeta[i] * ( matHyper->getS() - matHyper->getCt() * Gamma.segment<6>(6*i) );
                }
            }
            else
            {
                for(int i=0 ; i<nTau ; ++i)
                {
                    // S += vBeta[i] * ( matHyper->getSiso() - matHyper->getCtiso() * makeVoigt(Gamma.segment<6>(6*i)) );
                    S += vBeta[i] * ( matHyper->getSiso() - matHyper->getCtiso() * Gamma.segment<6>(6*i) );
                   // cout << "Siso= " << endl << matHyper->getSiso() << endl;
                   // cout << "Ctiso= " << endl << matHyper->getCtiso() << endl;
                }
            }
            break;
                
        case true:

        	for(int i=0 ; i<nTau ; ++i)
            {
                calcSrhe(i,C);
                calcCtrhe(i,C);
                S += Srhe - Ctrhe * Gamma.segment<6>(6*i);
            }
            
            break;
    }
}

// ***

void Viscoelastic::calcQv(const int & i, const VectorXd & C, const VectorXd & Gamma )
{
    switch (rheological)
    {
        case false:

        	matHyper->calcS(C); // If decoupled both Siso and Svol are computed

            if(!matHyper->hypertype.compare("coupled"))
            {
            	Qv = - ( 2.0 * vMu[i] * Gamma.segment<6>(6*i) - vBeta[i] * matHyper->getS() );
            }
            else
            {
            	//cout << "Siso= " << endl << matHyper->getSiso() << endl;
            	Qv = - ( 2.0 * vMu[i] * Gamma.segment<6>(6*i) - vBeta[i] * matHyper->getSiso() );
            }
            
            break;
            
        case true:
            
        	calcSrhe(i,C);
        	Qv = - ( 2.0 * vMu[i] * Gamma.segment<6>(6*i) - Srhe );

            break;
    }
}

// ***

void Viscoelastic::calcCt( const VectorXd & C, const VectorXd & Gamma)
{
    matHyper->calcCt(C);
    
    switch (rheological)
    {
        case false:
            
            if(!matHyper->hypertype.compare("coupled"))
            {
                // Ct_vis = [1 + sum_a (Beta_a)] * Ct - sum_a [ Beta_a * ( d Ct / dC) : Gamma_a]
                Ct = ( 1.0 + sumbeta ) * matHyper->getCt();
                for(int i=0 ; i<nTau ; ++i)
                {
                	calcDcPsiGamma(i,C,Gamma);
                    Ct -= D;
                }
            }
            else
            {
                // Ct_vis = Ct + sum_a (Beta_a) * Ct_iso^inf - sum_a [ Beta_a * ( d Ct_iso / dC) : Gamma_a]
                Ct = matHyper->getCt() + sumbeta * matHyper->getCtiso();
                for(int i=0 ; i<nTau ; ++i)
                {
                	calcDcPsiGamma(i,C,Gamma);
                    Ct -= D;
                }
            }
            break;
            
        case true: // Ct = Ct + sum_a (mu_a) * 1
                
            Ct = matHyper->getCt();
            for(int i=0 ; i<nTau ; ++i)
            {
            	calcCtrhe(i,C);
                calcDcPsiGamma(i,C,Gamma);
                Ct += Ctrhe - D;
            }
            break;
    }
}

// ******************************* RHEOLOGICAL  RESPONSE ********************************************

void Viscoelastic::calcWrhe(const int & i, const VectorXd & C)
{
    if(!matHyper->hypertype.compare("coupled"))
    {
    	 E = 0.5 * (C - Idv);
         Wrhe = vMu[i] * dotVoigt(E,E);
    }
    else
    {
    	matHyper->calcInvariants(C);
    	E = 0.5 * ( pow(matHyper->getJ(), -2.0/3.0) * C - Idv);
        Wrhe = vMu[i] * dotVoigt(E,E);
    }
}

//***

void Viscoelastic::calcSrhe( const int & i, const VectorXd & C)
{
    if(!matHyper->hypertype.compare("coupled"))
    {
        Srhe = vMu[i] * (C - Idv);
    }
    else
    {
    	matHyper->calcInvariants(C);
    	matHyper->calcP(C);
        VoigtFormat(matHyper->getP(),PMat);
        Srhe = vMu[i] * pow(matHyper->getJ(),-2.0/3.0) * PMat * (pow(matHyper->getJ(), -2.0/3.0) * C - Idv);
    }
}

// ***

void Viscoelastic::calcCtrhe(const int & i, const VectorXd & C)
{
    if(!matHyper->hypertype.compare("coupled"))
    {
        Ctrhe = vMu[i] * MatrixXd::Identity(6,6);
        Ctrhe(3,3) *= 0.5;
        Ctrhe(4,4) *= 0.5;
        Ctrhe(5,5) *= 0.5;
    }
    else
    {
    	matHyper->calcInvariants(C);
    	matHyper->calcP(C);
    	matHyper->calcPt(C);
    	matHyper->calctildeP(C);

    	// Jm = pow(matHyper->getJ(), -2.0/3.0);
        VoigtFormat(matHyper->getP(), PMat);
        VoigtFormat(matHyper->getPt(),PtMat);
        VoigtFormat(matHyper->gettildeP(), tildePMat);
        
       // cout << "Sa= " << endl << Sa << endl;
        
        calcSrhe(i,C);
        Ctrhe = 2.0 * vMu[i] * pow(matHyper->getJ(),4.0/3.0) * PMat * PtMat - 2.0/3.0 * (Srhe * makeCiVoigt(C).transpose() + makeCiVoigt(C) * Srhe.transpose()) + 2.0/3.0 * vMu[i] * Jm * dotVoigt(C,Jm * C - Idv) * tildePMat;
    }
    
    // cout << "Ctrhe= " << endl << Ctrhe << endl;
}

// ***

void Viscoelastic::calcDcPsiGamma(const int & i, const VectorXd & C, const VectorXd & Gamma)
{
    if(!matHyper->hypertype.compare("coupled"))
    {
        switch (rheological)
        {
            case false:
                D.setZero(6,6);
                break;
            case true:
                D.setZero(6,6);
                break;
        }
    }
    else
    {
        D.setZero(6,6);

    	matHyper->calcInvariants(C);
    	matHyper->calcCTensor(C);
    	matHyper->calcCiTensor(C);
    	matHyper->calcP(C);
    	matHyper->calcPt(C);
    	matHyper->calctildeP(C);

        VoigtFormat(matHyper->getP(), PMat);
        VoigtFormat(matHyper->getPt(), PtMat);
        VoigtFormat(matHyper->gettildeP(), tildePMat);
        InverseVoigtFormat(GamTens, Gamma.segment<6>(6*i));
       //cout << GamTensor(1,1) << "\t" << GamTensor(2,1) << "\t" << GamTensor(3,1) << endl << "\t"  << GamTensor(2,2) << "\t" << GamTensor(3,2) << endl << "\t" << "\t" << GamTensor(3,3) << endl;
       // cout << "Gamma= " << endl << Gamma.segment<6>(6*i) << endl;
        
        if(tildeGamoCi.size() == 0) tildeGamoCi.setZero(6,6);
        if(CiotildeGam.size() == 0) CiotildeGam.setZero(6,6);
        if(tildeGam.size() == 0) tildeGam.setZero(6);

        VoigtFormat(-1.0 * odot(matHyper->getCiTensor(), doubleContraction(selfOdot(matHyper->getCiTensor()), GamTens)), CiotildeGam);
        //cout << "CiotildeGam= " << endl << CiotildeGam << endl;
        VoigtFormat(-1.0 * odot(doubleContraction(selfOdot(matHyper->getCiTensor()), GamTens), matHyper->getCiTensor()), tildeGamoCi);
        //cout << "tildeGamoCi= " << endl << tildeGamoCi << endl;
                
        VoigtFormat(-1.0 * doubleContraction(selfOdot(matHyper->getCiTensor()), GamTens), tildeGam);
        //cout << "tildeGam= " << endl << tildeGam << endl;
        
        switch (rheological)
        {
            case false:
                
                break;
        
            case true:
                
            	calcSrhe(i,C);
            	calcCtrhe(i,C);
                D = - 4.0 / 3.0 * vMu[i] * pow(matHyper->getJ(),-4.0/3.0) * (2.0 * dotVoigt(makeCiVoigt(C), Gamma.segment<6>(6*i)) * PMat * PtMat + (tildeGam * C.transpose()) * PtMat + (makeCiVoigt(C) * Gamma.segment<6>(6*i).transpose()) * PtMat + PMat * (Gamma.segment<6>(6*i) * makeCiVoigt(C).transpose()) + PMat * (C * tildeGam.transpose()));
                D += 4.0 / 3.0 * vMu[i] * pow(matHyper->getJ(),-2.0/3.0) * (( 2.0 * dotVoigt(PtMat * Gamma.segment<6>(6*i), pow(matHyper->getJ(),-2.0/3.0) * C) - dotVoigt(PtMat * Gamma.segment<6>(6*i), Idv)) * tildePMat + dotVoigt(C, pow(matHyper->getJ(),-2.0/3.0) * (C - Idv)) * (tildeGamoCi + CiotildeGam - 1.0/3.0 * (tildeGam * makeCiVoigt(C).transpose() + makeCiVoigt(C) * tildeGam.transpose())));
                D -= 2.0 / 3.0 * ((Ctrhe * Gamma.segment<6>(6*i)) * makeCiVoigt(C).transpose() + 2.0 * Srhe * tildeGam.transpose() + 2.0 * tildeGam * Srhe.transpose() + makeCiVoigt(C) * (Ctrhe * Gamma.segment<6>(6*i)).transpose());

                break;
        }
    }
    //cout <<"D= " << endl << D << endl;
}

// ***

VectorXd Viscoelastic::getTau()
{
	VectorXd Tau(nTau);
	
	for(int i=0 ; i< nTau ; ++i) Tau(i) = vTau[i];
	
	return Tau; 
}

// ***

VectorXd Viscoelastic::getBeta()
{
	VectorXd Beta(nTau);
	
	for(int i=0 ; i< nTau ; ++i) Beta(i) = vBeta[i];
	
	return Beta; 
}

// ***
